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Techniques for extrapolating statistical data on the parameters 
of individual components to the statistical performance index for an 
overall system are considered. Two cases are SARIN (l) The 
deterministic case in which the system's performance index is known 
functionally in terms of the system parameters. (2) The non- 
deterministic case in which only limited data on the performance index 
and its sensitivity with respect to system parameters is known. Com- 
puter programs are developed in both cases for combining given 
probability density distributions for the parameters into an overall 
probability density distribution for the system performance index. 
Theory and programs are developed and verified with a specific numer- 


ical example. 
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lL. INGRODUCTION 


The performance ofa system, from a statistical point of view, 
depends upon the nominal values of the parameters of the component 
parts included in the system and on the statistical distributions of 
those parameters. Alternately, the response of the component parts 
to a given signal may be known statistically and it is desired to pre- 
dict an overall statistical distribution for the system response. 
Statistical information on the system's response may be interpreted 
in terms of probability of success or failure (operation or non- 
operation) for a single system or in terms of the percentage of 
operational systems that can be expected in an ensemble of systems. 
Currently three basic methods are being used to obtain performance 
data for systems; the worst-case method [Mark and Stember 1965], 
the moment method [Mark and Stember 1965], and the Monte Carlo 
method [Shreider 1966]. 

The worst-case method is a non-statistical approach intended to 
determine whether it is possible, with given parameter tolerance 
limits, for the systems performance characteristics to fall outside 
specifications. This method determines nothing of the probability 
distribution of the performance. Worst-case analysis, although it 
may be involved, does provide useful engineering insight and is used 


a great deal in practice. However, its outlook is basically limited. 


Statistics are combined with system analysis techniques in the 
moment method to estimate the probability that performance will re- 
main within specified limits. The actual parameter probability 
distributions are usually replaced by Gaussian distributions having 
the same variances. Hence, the moment method has inherent limita- 
tions. Parameter probability distribution characteristics such as 
skewness and peakedness are usually not taken into account. 

In the Monte Carlo method a large number of alternate replicas 
of a system are simulated by mathematical models. Component 
values are randomly selected and the performance of each replica is 
determined for its particular set of components. The accuracy of the 
Monte Carlo method is dependent upon taking a sufficiently large 
number of replicas and can be both lengthy and costly. Monte Carlo 
methods are also used extensively in practice although they provide 
limited engineering information. 

The purpose of this thesis is to present a numerical technique 
for combining statistical data on several independent components into 
a single distribution function for the system performance. This is 
accomplished by considering two general cases. The deterministic 
case, where the functional relationships between the system parameters 
and the output are known functionally; and the non-deterministic case, 
where these functional relationships are unknown but measured dis- 
crete data is available. In both cases the actual parameter prob- 


ability distributions are used in conjunction with numerical integration 
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techniques to obtain the system distribution function. One major 
advantage of the approach is that the known parameter distributions 
are quite arbitrary and may be histograms obtained from the results 


of experiments on individual components. 


li 


lige E Pe RIMINIS WIC ~CASE 


The object of the deterministic case is to investigate the dis- 
tribution of a system output when the output has a known mathematical 
relationship to several internal parameters. The basic assumptions 
are that all of the internal parameters are independent and that the 
probability distributions of the internal parameters, although arbi- 
trary, are known mathematically or given in terms of histograms. 

To commence the investigation, the general equations for various 
combinations of parameter distributions which may be expected ina 
typical system performance index must be obtained. These will | 
basically consist of operations on the distributions of random variables 
which are interrelated by addition, subtraction, multiplication, and 
division. The basic mathematics for these operations are well known 
in the literature, at least in general form [Beckmann 1967]. 

Consider two independent random variables: X, assuming values 
x with a known probability density py(x), and Y, assuming values y 
with a known probability density py(y). Also consider the functions 
U = U(X, Y) and V = V(X, Y) and their inverse functions X = X(U, V) 
and Y = Y(U. V). Then, when X is near x and Y is near y, U and V 


must be near u = u(x, y) and v = v(x, y). Hence, 


Pyy(*, y) dxdy = pyyy(u, v) dudv (1) 
Substituting the inverse functions and the required probability 


densities, we obtain 
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Resi 


where pyjy(u, v) is the joint probability density of U and V, Pylx(u, v)] 


Pry lu, v) = Bee [x (u, v)] Pe ala Ly (u, v) 


is the probability density of X, py [y(u, v)] is the probability density 





of Y, x(u, v) and y(u, v) are the inverse functions, and arn - 7 
? 


is the Jacobian, given by 


ax ak 
core _ |au ov (3) 
ee ay 2a 

au oV 


When X and Y are related by a single functional the foregoing can be 


modified by letting U = U(X, Y) and V = X. (4) 
Then, 
2X 1 
o( x _ | au _ . ay (5) 
(u,v ay A gu 
ou 


Equation (2) then becomes 


Puy (4, oo Py (x)Py lec. 22 ie (6) 
co 
Since, in general, Py(a) = | Pyx (u,x)dx, equation (6) may be written 
= 00 
as 
66 
Pylu) = J pxvery [ve x)] |22| ax (7) 
0 


Using equation (7), the specific formulas of Table I have been 
compiled for the functional relationships and operations listed, which 
were chosen on the basis of what might be expected in engineering 


problems. 
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The technique outlined here consists of the successive application 
of the formulas of Table las required, with the required integration 
performed numerically. Various types of numerical integration are 
possible; rectangular integration, trapezoidal integration, or Simpson's 
Rule. These are summarized as follows: 


Rectangular integration [Stanton 1961], 


1s ri 
Jets i = AF fia) (8) 


Xx 7 =: 


Trapezoidal integration [Stiefel 1963], 


Xn Ned 
Oa | a 
fecnex AX ar (x, at (x,,)| + £(x, ) (9) 
a2 


Simpson's Rule [Stanton 1961], 


X,+naXx tL 
_— As: 
ie R/Gz =: ae aS, ) oe 


Xo IO 


where C = 1, 4& & --- 2, 4, 1 

Although the specific technique of numerical integration is of 
prime importance in so far as accuracy and computer programing 
and time are concerned, the exact procedure used is a detail which 
is best dictated by the particular problem to which these techniques 
are to be applied. 

This method may be best illustrated by considering the following 
example. Consider that: 


H =  - b 4 Z (11) 





where a and b are constants and X, Y, and Z are independent random 
variables with known (measured) probability densities. Consider that 


Py(x), Pyty); and P7(z) are given as follows: 


meet x2) Py (y) ae) 
2 2 2 
abe lk 
O65 O45 1 O45 1 
X y Z, 


I 
) 
-~ 


(a) Consider first the combination U = XY, with gx 


Using rectangular integration for simplicity, equation (3) of Table I 


becomes, for w = Uy 
20 - (x.) Te 
pyla,) = [py(x)py(u/x) 8 = oxy en] (aliZ ) 
= =i * 
When u= u, = 0.09 
By (x: A 
Xi Ui /x; 2) R (41/x;} fx) R (4:/y:) 
0.1 0.9 4.0 1.8 auc 
O22 0.45 4.0 0.9 3.6 
OF 3 0.3 4.0 OY 6 2.4 
0.4 Graz 2 4.0 0.43 lL. we 
@. 5 0.18 4,0 O37 I. 43 
0.6 0.15 2256 0. 30 0. 78 
Oy? Oeiks ee 0.26 0. 44 
0.8 0.11 ] #0 Ov. 21 O21 
0.9 0.1 0.44 0.20 0. 088 





sum= 17.918 
Py (Oni l7- 918) = 1.79 


bs 


Be(%:) R(X) 
<-_* BRIDES U 
Xi U2 /x; Xi Ry (U2/x; ) Xs pi 2/y;) 
Diz) Zz. 0 4.0 0.0 0.0 
Ore2 ew 4.0 Z20 8. 0 
OS Oi. 4.0 1.4 ore 
0.4 Sead: 4.0 0 4.0 
Ones: 0.4 4.0 0.8 36 
0.6 Ors3 Zc O26: / 1.74 
On OS 1.7 Oar P02 
0.8 ORZ 5 0 0e5 0555 
On S, 0.22 0.44 0.42 Oe les 
sum = 24.245 
Pylu,) =mO. 1) (24.245) = 2.4 
The pattern for evaluating p,,(u) is repeated for u=.3, .4, --, 1.0 
yielding the curve of Figure 1. 
(b) Next consider the function V = i = im - Here equation 
U XY 
(7) of Table I may be applied. 
“SPy, (iv) = go Cn } ( 3) 
x te ) ’ 
; = - | A _ J. ee 
Deed 10.0 eo 0.01 Ol Olas 
Olea a0 Zone 0. 04 Bhs 0% 
O43 = 33 2259 O08 Oi215 
0.4 ae5 2.8 0.16 0. 338 
05 Zid 1.54 O25 0.385 
0.6 Il aaex henO2 0 0. 364 
Opt 1, 43 Je55 0.49 02270 
0.8 Nerage: 0.28 0.64 ORT 
On? re O05 Cell 0. 040 


Py(v) is plotted in Figure 2. 
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(c) The function W may now be calculated using equation (8) 





of Table Ir where W = aVtb= = oi (14) 
and 
5 ee ~~ “ae (= a a “al by ( a y ( l 5) 
If a = 2andb = 1, the following result is obtained. 
V W=aVtb Py(v) Py(v)/2 
l 5 0.0 0.0 
és 5 O7S55 071925 
3 7 G2 25 ©. #25 
4 9 0.15 Oe 075 
5 ie 0.10 0.050 
6 ik 0. 06 0.030 
% 15 0. 04 0. 020 
8 17 0.03 0.015 
9 19 0.02 0.010 
10 Fil 0.01 0.005 


The result is plotted in Figure 3. 


(d) Since, from equations (11) and (14), 
H= Wt Z (16) 
it follows from equation (1) of Table I that p;,(h) is given by the con- 
volution of py(w) and Py (z). This operation may also be carried out 


numerically and for the given p7(z) may be expressed as 


h 


pi(h) « |a-(e)de for = be > 
EY ‘ 
0 
(17) 
hy 
: olds) = fey ws for n> i 
h-| 


These integrals may also be carried out numerically yielding the 


results plotted in Figure 4. 
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To work more complicated problems than the previous example 


requires a great amount of mathematical calculations. 


Thus, when 


using numerical integration techniques the method lends itself nicely 


to the digital computer. 


graphs corresponding to some of the equations of Table I. 


Figures 5 through 9 are the computer flow 


The com- 


puter program for any mathematical function can then be derived 


from the appropriate combinations of these basic flow graphs. 


To illustrate a more complex problem consider the following 


example: 
given 
© = 1285 ES - = + W 18 
A B (18) 
where 
A = nan ara B = B3 Bey 
and 
a B3pb 
3a 1 + Bay(2. 16 x 10-2) 
— P 6b 
6a = 
b + Be, (0. 63% 10-4) 
Bay, is a normal distribution; mean = 130, standard deviation = 10 
Bg, is a normal distribution; mean = 65, standard deviation = 12 
W is a uniform distrmsution from te- - 4/7 to wet 4/7 


Theta is system output parameter and the variables B3,, Bop: and 


W are independent. 
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Since the variables are independent, equation (18) can first be 


simplified by ordinary mathematical manipulation. 


i 
8 


om is 2 % -2 
Let *& Y = Bey, ky = 2. oa TOS k> =" .162 ax FO 


S30) 


then B = XY 


- - == kik, t q as *2 
A B Y x 
and @ = 0.175 + 27.8 + 8.1 4 W (19) 
xy Xx 


Now, utilizing the basic methods and flow graphs, the problem 
can be programed on the digital computer and the probability density 


of 8 readily obtained. The program and solution are contained in 


Appendix I. 
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Ill. NON-DETERMINISTIC CASE 


In many practical problems the functional relationship between 
variables is unknown and it cannot be expressed simply by formula. 
It is generally available only as numerical data from a computer run 
on a mathematical model of the system under study. Two cases are 
considered here. First, the case where the system response of 
interest, U(X), and its sensitivity coefficients, U'(X), are available 


at discrete values of a single variable, X. 


U(x) = [uloy), ule), - - -, ulx,)] (20) 
U(x) = [ul(x,), ule), - - -, u(x,)| (21) 
Second, the case where the system response of interest, U(X. Y), 
and its sensitivity coefficients, U'(X, Y), are available at discrete 


values of two independent variables, X and Y. 


palgleredia 6 6 iy 
Dix, ¥) = U2] nD Sa oer 
(22) 
el Cae = Ue 
where 
1j a u(x. , a 
and 


20 


Uingex., Y) = u'21 u ee - u 
(23) 


where 


In both cases the problem reduces to determining p,,;(u) from the 
known distributions Px (x), Pyly) and the data of equations (20) and (21), 
Or (22) and (23). 

In the first case, consider that the probability density, Py (x), 


and the cumulative distribution, 
Ky 
ak 5 


- 0 


are known. Then at each data point of u(x), 
Re | ule;) = SP ileg) (25) 


That is, the probability that u <= eh is equal the probability that 
x = x.. Equation (25) gives discrete points on the cumulative 
distribution, Pi;(u), for each value of the parameter x for which 
u(x) has been measured. 


It is also well known [Beckman 1967] that 


"~" fs, ‘ 
A: fy \ se? 


he. ae | =) (26) 


ta 


Thus the slope of the cumulative distribution curve, -- 50," lm 


at each data point is equal to the known probability density, Bees 


cu 
divided by the measured sensitivity function, |- ~ , at the 


7” 
- 








& 


data point. The conclusion is that the cumulative distribution, Pilu), 
and its slope may be calculated readily from the measured response 
and sensitivity data. This curve is therefore evaluated as a piece- 
wise linear approximation. 

When U is a function of two variables, X and Y, the foregoing 
cumulative distribution curve can be evaluated for each value of the 
second variable, Y. Thus a family of such curves or data points 
denoted by P_,(u, y;) is available. For a specific value of u = u; 
each of these curves gives the probability that u =u, when y = Yj: 
The total probability that uS u, depends uponpy(y)Ay and is 
given by the equation: 


m 
Py,(u;) = ay) Pala, Vs) Py(v,) (27) 


e! 
foo 


where P,,(u, y;) is the cumulative distribution evaluated with y asa 
constant and x as a single variable as discussed previously. 
This result can be proven by considering the following (see 


equation (7) ): 


co 
BR 
2 


Py(u) = yi¥) 2x [2(u,y)] | Fey (28) 


-* 


where u = u(x, y). ow is the sensitivity function with respect to 
the parameter x and in equation (28) must be expressed in terms of 


uandy. It is convenient to denote 


ZZ 


Ca (u,v) 
a Gel (29) 


and then equation (28) becomes 


0% 
(a) = poy ly )eg| x(usy)} o-fe 30 
Py io By li i Py ALU, a (u,v) ( ) 
fe 


= 6 
We also know 
U 
Py (u) = fry(oiae (31) 
= 06 


Substituting equation (30) into equation (31) and inter-changing the 


order of integration yields 


U 


~ 
Puy (u) = fe: (yv) ssf, [> ( Us¥ ) a GI (32) 
J Us 4 
Now 

u u 

| Px [eos gees = fs (ao dédiarsP ppp laaeat) (33) 
Age ay ) U ——— 

— 0 tise 


Since, if Y is considered a constant and X as the random variable, 


u = u(x, y), 


_ Px lx yy] py [x(x y)] (34) 
pylu, y) = du uy(u, y) 
ox 


Substituting equation (33) into equation (32) yields: 


co 


Piya) = feyoreyt, y)dy (35) 


— 0 


In discrete form, equation (35) becomes equation (27) when 
rectangular integration is used. However, more accurate integration 


techniques may also be applied to equation (35), 


oo 


To demonstrate the application of this non-deterministic method 
consider part of the previous example of equation (18). It is now 
assumed that the functional relationship between the output and the 
parameters is unknown. However, we do know the probability den- 
sities of the parameters x and y and have certain measured data which 
consists of measured values of x, y, and § and the corresponding 
sensitivities for each measured %. The data is presented in 
Table II. Following the procedure outlined, the problem can be 
programed on the digital computer and the probability density of # 
readily obtained. The program and solution are contained in 
Appendix II. The results can be compared with those obtained using 
the deterministic formulas. Noting, of course, that the solution 
to the deterministic problem,contained in the Appendix, includes 
the addition of the uniform distribution W(w) and the solution to 


the non-deterministic problem does not. 
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IV. CONCLUSIONS 


The numerical technique presented here has distinct advantages 
over the worst-case, the moment, and the Monte Carlo methods. By 
utilizing the actual parameter probability distributions the results 
obtained are more accurate than those of the moment method and pre- 
sent a great deal more information than the worst-case method. The 
accuracy obtained is equivalent to that of a large number of Monte 
Carlo replicas. Thus the saving in time over the Monte Carlo method 
is obvious. In general, the numerical method presented here can 
obtain a more accurate solution in less time than any of the other 
current methods. However, certain limitations are inherent in the 
basic assumptions used to formulate the method. Since independent 
parameters are required, a great number of practical problems can 
not be handled by this method as it is presented here. By extending 
the mathematical derivations to include the correlation of parameters, 
the method could be used for systems in which the parameters are 


not independent. 


Le 


Table I 


Functional Formulas 


Ce a oe api) = { Py (x) py (u-x) ax 


= J Px (ary) py (yay 


(2) Uv= X so¥ Py(u) = f Px (x) py (x-u) ax 

= f py(uty)py(y)dy 
(3) U = XY pyla) = {Py (x) py (w/a) 
A Me Cee Xe Py(u) = 5 [Pg (x) Py (x/u) xt ax 

UP igh 
= { Py (uy )py(y) ly ay 
(a (V8 
(HyHU"'= x? Py (a) _ BS Bee: u20 
2vu 

C6) Us| HI Py(v) = Py (x) + Py (-x) > 0 
(7) U <—.a/k Py(u) = “2 Py (a/u) 


Ge) le) 


be 


1 - 
by = gr py(2 


(a and b are constants) 
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Table II 


Sample Data 
x J 41 oo 65 TT 89 
120 0.7669 0.6699 0.5554 


) 
¢ 


[hh «4 53 6 | 77 #| 489 


fe 





0.0165 | 0.00988 | 0.00658 { 0.00469 | 0.00352 


X 110 120 130 140 150 


ha 0.000670 | 0.000561] 0.000478 | 0.000413 | 0.000359 
Xx ) 


on 





Figure 1 





re) 8 FF TIS Ger 6s of 110811 
V 
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10 


Figure 3 
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Figure 4 
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1 


W= 2V +1 

5 20 

H = W+ 2 
20 


15 





Calculation of u, and p(u,) where U = a/X and a, x, 
9 


and p(x,) are known. 


Figure 5 


350 


START 






Read 
x5, p(x,), 4x 


i=l »N 















Compute 
p(x, by 
Interpolation 





al =. 


To perform rectangular integration of the distribution p(x, ) 
VB. X;; the distribution must be approximated by rectangular 
areas of equal widths ax. This program calculates the x's 
such that ve! = x. + Ax and interpolates for the correspond- 
ing values p(x.) from the given distribution p(x, ) v8. X,. 


Figure 6 


ps 


START U =X +Y 






Read 
Kir P(xy), ys Ply,) 
Ox -(i=1,N;j=1.M 


u(1)=x(1)+y( 


u(i)su(i-l)+ax 


Yes 


Coy a 


No 





a 
p(u;)=p(u,)+p(x,)*P(y 5443 )*o% 





This program performs the convolution U = X + Y by using 
graphical convolution techniques. The distributions p(x, ) 
VS. X,; and p(y) vs. y. are represented by sets of rec- 
tangular areas of width ax. The distribution of xX, is 
reversed and passed across the y. distribution in Ax in- 
crements. The distribution p(u, vs. u, is then obtained 
by the normal graphical convolution mathematical calcula-~ 
tions[Bracewell 1965, Healy 1969]. 


Figure 7 


29 





A cs. 





START 


Read 
px(xq), [*i]p , AX 
Py(yk), Yk» Uj 


Compute 





En y 








Compute 
Py(y,) by 
Interpolation 
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Py (y})=py(y;,) 


Form 
Py(y})| 





Compute 


R= fepPy (ys 








Yes 









Compute 


T = py(x,)R 





Figure 8 
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a an -- 






Compute 
y= = x5 /v, 


U = X/Y 










Py (y})=Py(y,) 








Compute 


Py(y}) by 
Interpolation 
Write 
Py(y}) 
Form 
Yes No 
| gaze 
p (v3) 
mS! Z 


Compute Compute 


S = axu;* X= Py (x) [xp 











Compute 


T =,K, py(yi) 





Figure 9 
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APB NDIx I! 


The solution to the deterministic problem has been obtained by 
using three separate programs. Linear interpolation and rectangular 
approximations have been used to perform all convolution and rec- 
tangular integration. 

The problem and given information is as follows: 


Obtain the cumulative distribution of 8 where 


a le 7 ay me $ (rw) 
ven Y(y) ' 

C = 0.175 

A = 8.1 

B = 27.8 


pix) ¥s IN( 130; 10) 
Y(y) is N(65;12) 


W(w) is uniform 0. 875 from - 4/7 to + 4/7 


Program one inverts the distribution x vs. p(x) and multiplies 
it by the constant A. Program two does the same for B/Y(y). Since 
both X(x) and Y(y) are normal distributions, the input data for these 
programs was taken from standard normal tables. 25 data points 
were used in each, the median and each 1/4 standard deviation out to 
3 standard deviations. The resulting distributions, yetS) and Y'(y), 
are then used as input data for program three. Program three per- 


forms the convolution U(u) = X (x) + Y'(y), the addition U'(u) = U(u) + C, 
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and the convolution 9 = U'(u) + W(w) to obtain the cumulative dis- 
tribution of 0. 

To obtain highly accurate results a Ax of 0.002 was used for 
convolution and rectangular integration. The over-all accuracy for 
this type of solution will depend upon how closely the rectangular 
approximations resemble the continuous distributions which they 


represent. 


ay, 
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The tabulated sample data for the non-deterministic case, Table II, 
was obtained from a modified form of equation (19). Using the 


modified equation, 





where X(x) is N(130;10) 


Y(y) is N(65;12) 


values of @ and sensitivity were calculated for chosen values of x and 
y. Assuming these to be measured values, the data was then used to 
obtain the cumulative distribution of @. 

Since the data indicates g is more sensitive to changes in y than 
to changes in x, the problem is worked with x as a constant and y as 
a single variable as discussed in section III. Aitkin-Hermite inter- 
polation [Hildebrand 1956] is used to estimate the values of y between 
data points. This interpolation scheme utilizes the sensitivity data 
and obtains a more accurate estimate than would linear interpolation. 
‘The caleulation of the cumulative distribution of y is then made 
using equation (35) in discrete form. As inthe deterministic case, 


a Ax of 0.002 was used for rectangular integration. 
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